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airfoil  and  the  results  were  compared  with  the  available 
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transition  model  has  an  important  influence  on  the 
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I.   INTRODUCTION 

There  are  many  theories  to  aid  the  airfoil  design 
process  by  computational  methods,  because  of  the  desire  to 
reduce  the  number  and  cost  of  wind  tunnel  tests.   A  still 
largely  unresolved  question  is  the  problem  of  flow  sepa- 
ration.  Because  the  classical  boundary  layer  approximation 
cannot  be  applied  to  separated  flow  calculations,  engineers 
have  tried  to  overcome  this  limitation  by  developing 
viscous/ inviscid  interaction  approaches  or  to  develop  direct 
solutions  of  the  Navier-Stok.es  equation. 

The  purpose  of  this  thesis  is  to  demonstrate  the 
capability  of  the  viscous/ inviscid  interaction  method  by 
applying  Cebeci's  interactive  computer  program  to  a  single 
airfoil  ( FX  63-137)  at  three  low  Reynolds  numbers  and  by 
comparing  the  results  with  experimental  data. 

Chapter  II  explains  the  boundary  layer  theory.   The 
boundary  layer  equations  are  derived  and  the  turbulence 
models  are  introduced.   Also,  this  chapter  includes  the 
prediction  of  transition  boundary  layer  calculations  and 
flow  separation. 

Chapter  III  introduces  the  interaction  methods. 
Three  weak  interaction  methods  are  explained  briefly  and  the 
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simultaneous  method  is  presented  as  a  strong  interaction 
method. 

Chapter  IV  describes  Cebeci's  interactive  computer 
program.   Input/Output  data  description  and  JCL  files  are 
included.   Also,  the  results  of  the  application  of  Cebeci's 
program  are  discussed. 
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II.   BOUNDARY  LAYER  THEORY 

A.   DERIVATION  OF  BOUNDARY  LAYER  EQUATIONS 

Generally,  the  thickness  of  the  boundary  layer  increases 
with  viscosity,  or  it  is  possible  to  state  that  it  decreases 
with  viscosity,  or  it  is  possible  to  state  that  it  decreases 
as  the  Reynolds  number  increases.   From  exact  solutions  of 
the  Navier-Stok.es  equations,  it  was  seen  that  the 
boundary-layer  thickness  [a)     is  proportional  to  the  square 
root  of  kinematic  viscosity  (v^); 


/ 


( X  )  «^ 


U 

where  x  is  the  distance  from  the  leading  edge  of  a  flat 
plate.   Using  the  local  Reynolds  number 


Re  =  U  x/i;,      /(x) 


(.-o 


For  simplicity,  assume  a  two-dimensional,  steady 
constant  -  property  flow  without  body  forces  and  leave  the 
stresses  unspecified  so  that  the  results  apply  to  laminar  or 
turbulent  flow.   If  simplifications  are  to  be  introduced 
into  Navier-Stokes  equations,  it  must  be  assumed  that  the 
boundary  layer  thickness  is  very  small  compared  with  a 
representative  linear  dimension  (L)  of  the  body,  ie.  cT  < <  L. 
In  this  way,  the  solutions  obtained  from  the  boundary  layer 
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equations  are  asymptotic  and  apply  to  very  large  Reynolds 
numbers . 

From  the  Navier-Stok.es  equations, 

x-direction:  -  1_  d  P  +  1^  d^xx    +  1^  ^^Icy  =  u  d  u  +  v  ^  u 

/*'ax/'ax   /'ay  ^x     ay 

y-direction:    -    1_     d  P    +    1^     ^6\v    +    1_     3<^xv    =   u  ^  v  +    v     d  v 

/"^y/^ay   /^ax  ^x     ay 


and  the  continuity  equation  is 
j>  u  +  a  V   =  0 

ax        ay 


(2.3) 


Inserting    the    "typical"    (order-of-magnitude )    values,    replace 
the    dependent    variables    as    follows; 


u    «^o    u. 


d  u  ue 

■  c-> 


a  u 


c^ 


ue 


ay        /  ax       '        i 

Then  Eq .  (2.1)  can  be  expressed  as 

-  i_   a  p  +  1  a<rxx  +  i_    ^(Txv  =  u  d  u  +  v  a  u 
/^ax^ax   /«ay      ^x     ay 


11  11 


ue 


(2.4) 


1 


where  the  typical  values  were  written  below  the  terms  to 
which  they  correspond. 

Considering  turbulent  flow,  all    stresses  are  of  the 
same  order  ie.  a  general  stress    must  be  of  order 
P\ie(f/l,    then  Eq.  (2.1)  becomes 


-  1  a  P  +  1  dtTxv 

/^  ax  p     ay 


1   +  Q(  <f 


if) 


=    u  du+v  au 


(2.5) 


a  X 


ay 
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where  0(<7/l)  indicates  a  quantity  of  the  order  of  magnitude 
o  f  //  1  . 

Considering  a  laminar  flow  of  Newtonian  viscous 
fluid. 


6"  XX  =  lyU   d  u   (Txy  -Xil  d  u  +  3  v^  (Tyy  =  2//  by 
^  X,          [ay    3x  /  ^y 

In  the  (Txy  term,  Ti^lb^i    is  smaller  than  3u/ay. 
Therefore,  Eq .  (2.5)  becomes 


-  I  _i_P  + 

/^  a  X 


V 


d^U 


a/*  L 


+  ofcf  ]" 


( —)       ~       ^      d  u  +  V  d   u 
U  ^  J       ax      dy 


Similarly,    Eq .     (2.2)    can   be    written   as: 

-    1^     3  P    +    1_      a  (Tyy    +    1^      3  (Txy    =    u      3  V    +    V     b  V 

/^dy/^dy       /^ax  dx  ay 


(2.6) 


(2.7) 


1     ^  P  Pue  /   Ue      Ue/\      u//  

/<>    dy  1/     ^l  1/   '    /W       1"  1" 


ue  d 


(2.8) 


(2.9) 


If    we    write    all    the    viscous    terms    together, 

\_     _d_P  UgV       jJ_       I      if         lrl£_f  U//  u// 

/^     a  y         1-"       uel     I  //    L      'l    '    J  1""  l"" 

It    is    known    that    [cf/D   t^l^/n^l    is    laminar    flow,    so    that    the 

J-      y 

viscous  terms  are  also  of  order  Ugd/l   ie.  bP/^y    is  of 
order  u///l  ,  but  the  pressure  difference  between  y  =  0  and 
y=  d    is  of  order  pag  d/1    and  the  difference  in  5P/dx  will  be 
negligible  compared  to  the  external  stream  dynamic  pressure, 
1/2  pMg    .       Therefore,  for  practical  purposes. 


i_P  =  0 

ay 


(2. 10) 
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For  this  case,  since  changes  in  P  must  be  of  the  same  order 
as  changes  in  y,  the  pressure  does  not  change  significantly 
through  the  boundary  layer. 

Thus,  the  entire  equation  of  motion  in  the 
y-direction  may  be  dropped  from  further  considerations.   In 
this  way,  the  following  simplified  equations  are  left  for 
the  analysis  of  a  boundary  layer: 

_a_u   +  _a_v   =0  (2.3) 

■2. 

-    1_  d  p  +  ^  d   u   =   u  a  u   +   V  ^  u  (2.11) 

/^     d  X  d  y"^  d  X  ay 

a  p  =0  (2. 10) 


ay 

These  relations  are  known  as  Prandtl's  boundary  layer 
equations.   Unless  one  encounters  very  high  Mach  numbers, 
the  above  orders  of  magnitude  are  not  changed  when 
compressibility  effects  are  considered. 

B.   LAMINAR  AND  TURBULENT  BOUNDARY  LAYER 

The  low  viscosity  fluid  flow  past  solid  bodies  should  be 
considered  as  two  regions.   One  is  the  thin  region  near  the 
boundary  in  which  the  effects  of  viscosity  are  concentrated 
and  the  other  is  the  region  further  away  from  the  wall  in 
which  the  influence  of  viscosity  is  so  small  that  it  can  be 
neglected.   Thus,  it  can  be  stated  that  all  viscous  effects 
are  concentrated  in  the  thin  region  which  is  known  as  the 
boundary  layer.   This  boundary-layer  type  behavior  requires 
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high  Reynolds  numbers.   Generally,  the  thickness  of  the 
boundary  layer  (/)  is  defined  as  that  distance  from  the  wall 
where  the  velocity  (u)  differs  by  1%  from  that  ( Ue ) 
calculated  by  the  ideal  flow  analysis. 

Consider  a  constant-property,  steady,  two-dimensional 
flow  past  a  flat  plate.   If  u/Ue  is  plotted  against  a 
dimensionless  y-coordinate,  T  =  (Ue/vx)  V   ,  the  velocity 
profiles  are  geometrically  similar  and  reduce  to  a  single 
curve  for  a  laminar  boundary  layer  flow.   This  is  well  known 
as  the  Blasius  profile.   The  geometrical  similarity  is 
maintained  regardless  of  the  Reynolds  number  of  the  flow  or 
of  the  local  skin  friction.   For  a  turbulent  boundary  layer 
flow,  since  the  viscous-dependent  part  and  the  remaining 
Reynolds-stress-dependent  part  of  the  profile  require 
different  length  scaling  parameters,  there  is  no  choice  of 
dimensionless  y-coordinate  that  leads  to  the  collapse  of  the 
complete  set  of  velocity  profiles  into  a  single  curve. 

The  conspicuous  difference  in  profile  shape  between 
laminar  and  turbulent  shear  layers  can  be  found  in  the  wall 
flows.   Because  of  the  constraint  on  eddy  size  in  wall 
flows,  the  efficiency  of  turbulent  momentum  transfer 
decreases  rapidly  near  the  wall.   But,  the  efficiency  of 
viscous  momentum  transfer  is  not  dependent  on  y  distance  in 
the  flow  which  has  no  heat  transfer. 
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C.   TURBULENCE  MODELS 

The  unsteady  continuity  and  Navier-Stokes  Equations  are 
valid  in  both  laminar  and  turbulent  flows.   But,  it  is  too 
difficult  to  deal  with  the  instantaneous  properties  in 
turbulent  flow  because  the  turbulent  flow  has  a  complex 
time-dependent  behavior.   Thus,  the  following  turbulence 
models  are  used  to  make  the  analysis  of  turbulent  flow  more 
convenient. 

1 .   Prandtl's  Mixing-Length 

Consider  two  adjacent  stream  layers  of  fluid  which 
move  with  different  velocities.   If  a  particle  of  fluid 
moves  from  one  layer  to  the  other,  a  momentum  change  occurs 
between  two  layers.   The  fast  particles  which  enter  the  slow 
layer  make  it  faster  and  the  slow  particles  which  enter  the 
faster  layer  impose  a  drag  on  it. 

The  mean  velocity  of  a  stream  layer  is  u,  and  that 
of  the  other  is  u  +  1  d   u/^y  where  1  is  the  distance  between 
two  layers.   Also,  the  fluctuating  velocity  in  the 
x-direction  is  u',  and  that  in  the  y-direction  is  v'. 
Prandtl  assumed  that  the  turbulent  fluctuations  are  due  to 
the  difference  in  the  mean  velocities  of  the  two  layers.   So 
u'  =  1  au/ay  ie.  the  fluctuating  velocity  in  the  x-direction 
is  of  the  order  of  the  difference  in  the  mean  velocities  of 
two  layers  which  have  a  distance  1,  where  1  is  the  mixing 
length.   Prandtl  also  assumed  that  all  components  of 


17 


fluctuating  velocity  at  a  given  point  are  of  the  same  order 
of  magnitude.   Thus,  v'  can  be  defined  as  v'  =  kl  au/ay 
where  k  is  a  constant. 

The  turbulent  shear  stress  due  to  momentum  exchange 
between  two  layers  is  the  rate  of  momentum  transfer  per  unit 
area.   Then  the  mean  turbulent  shear  stress  on  the  fluid  is 
T^=  -^u'v'  where  u'v'  ^**^|l  au/^y|  |  kl  ^u/ay|  .   Since  the 
values  of  1  and  k  are  unknown,  combine  these  two  unknowns. 
Then,  Lf=^l  pu/aypu/ay  where  1  is  called  the  mixing 
length . 

2.   Cebeci-Smith  Model 

With  the  eddy  viscosity  concept,  the  momentum 
equation  for  2 -dimens ional  laminar  and  turbulent  boundary 
layers  can  be  written  as: 


{  b  f  "  )  •  +  m  +  1 


ff"  +  m[l  -  (f'fj   =  xff'_^'  -  f  "_^|.  \ 

^    dx        d  x  J 


where  x  is  the  transformed  x-variable 

i"  =  ^m/  V 

Id  =      (1  +  t)   (l+£) 

Let  the  turbulent  boundary  layer  be  a  composite 
layer  consisting  of  inner  and  outer  regions.  Then,  the 
eddy-viscosity  formula  for  the  inner  region  is: 


(  ^m  );  = 


'^  (t)' 


<>u 


ay 


rt 


(0  i  y  i  y^.  ) 


where  1  is  the  mixing  length  for  2-dimensional  flow 
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/. 


'   =   1  -  exp  -     __1      Rx,^  / 
^  L   1200    ^'^    ^ 


tr 


is  an  intermi ttency  factor  for  a  flat  plate. 
R^    =  Up    Xff/i)   is  the  transition  Reynolds  number. 
The  eddy-viscosity  formula  for  the  outer  region  is: 

fCO 


(^'T,).  =  -< 


( Uf  -  u )  dy 


ftr    (y  ^  y  i^) 


where  R^  >_  5000,  the  universal  constant  o<    =  0.0168 

R^  <  5000,  c^    varies  with  R^  according  to  the 
the  empirical  formula 


'/j 


oC     =  0.0168    1.55 
1  +  A 

A  =  0.55  [l  -  exp  (0.243  h  ^-   0.298K)] 

K  =  Re  /4  25  -  1 

But  this  model  is  not  used  in  Cebeci's  interactive  computer 

program  which  will  be  presented  in  Chapter  IV. 


D.   TRANSITION 

The  boundary  layer  with  a  finite  thickness  starts  out  as 
laminar  at  first  in  the  flow  past  an  airfoil.   However,  the 
boundary  layer  becomes  unstable  and  all  small  disturbances 
begin  transition  to  the  erratically  unsteady  condition  which 
is  known  as  turbulence. 

In  the  boundary  layer  on  blunt  bodies,  transition  makes 
the  point  of  separation  move  downstream  which  decreases  the 
width  of  the  wake.   There  is  an  abrupt  change  in  the  drag 
curve  of  a  sphere.   This  change  is  due  to  a  boundary  layer 
effect  and  is  also  one  of  the  transition  phenomena. 
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The  process  of  transition  on  a  flat  plate  at  zero 
incidence  shows  a  sudden  increase  in  the  boundary  layer 
thickness.   Furthermore,  transition  involves  a  great  change 
in  the  shape  of  the  velocity  distribution  and  a  large 
decrease  in  the  ratio  of  the  displacement  thickness  to  the 
momentum  thickness.   Also,  it  causes  a  large  change  in  the 
skin  friction. 

As  described  in  the  above,  the  transition  from  laminar 
to  turbulent  flow  plays  a  very  important  role.   Since  the 
transition  is  not  an  instantaneous  process  and  the  flow  is 
intermittently  laminar  and  turbulent  over  a  certain  length 
of  the  airfoil,  a  certain  point  where  transition  takes  place 
cannot  be  assigned.   At  present,  there  are  no  exact  methods 
to  calculate  the  flow  within  the  transitional  region. 
Nowadays,  however,  two  methods  -  Michel's  method  and  the  e 
method  -  are  used  to  predict  the  transition  empirically. 

1.   Michel 's  Method 

Michel  investigated  many  kinds  of  data  for 
incompressible  flows  without  heat  transfer  and  found  this 
empirical  correlation. 


Rq  =  1. 174 


1  +  22400 


R 


0.4-6 
fr 


+r 


where  R^  =  Ue  e/P,   R;j  =  Ue  x/p  . 

Since  the  momentum  thickness  grows  more  rapidly  when 
the  pressure  gradient  is  positive,  Michel's  equation 
involves  the  effect  of  pressure  gradient.   However,  the 
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effect  of  surface  roughness,  which  is  also  very  important, 
is  not  included. 

2.   The  e^  Method 

f 

This  e   method  is  of  spatial  amplification  theory 

based  on  the  solution  of  the  Orr-Sommer feld  equation  which 
forms  the  point  of  departure  for  the  stability  theory  of 
laminar  flows . 


g' 


=  a  =  esp 


^i   -  /^i  d  o^r 


) 


dX 


where  X  =  x/  ^ 


o('i  -  ^  ' 


I  oC\    (f 


d  /Sr 
^\  =  Idi     (f 


g'(x,y,z,t)  is  a  typical  disturbance  function. 

<:/-  -    2  J     (  A  is  the  wavelength  of  the  disturbance)  . 
A 
/^ r  is  the  circular  frequency  of  the  partial  oscillation, 

/^r  is  the  amplification  factor  which  determines  the 
degree  of  amplification 

The  basic  assumption  is  that  transition  begins  at 

the  point  where  a  small  disturbance  introduced  at  the 

critical  Reynolds  number  is  amplified  by  a  factor  of  e 


LAMINAR  BOUNDARY  LAYER  CALCULATIONS 

1 .   Similar  Solutions 

a.   Blasius  Solution  For  a  Flat  Plate 

Assume  -  a  flat  plate  at  zero  angle  of  attack 

-  steady,  2-D,  laminar,  incompressible 
flow . 
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-  constant  viscosity 

-  negligible  body  force 

Since  the  pressure  along  the  plate  is  constant, 
there  is  no  pressure  gradient.   The  simplified  Navier -Stokes 
Equation  is: 


u   3  u 

ax 


V  a  u  =  ,)  LH 


ay 


^J7^ 


(2.12) 


(B,C)    y=0;   u=0,   v=0 
y  =  CO  ;   u  =  Uoo 

To  transform  from  the  partial  differential 

equation  to  an  ordinary  differential  equation,  define  the 

following  transformation  parameter: 

1    =  /  Ua>    y        where    ^  =  f(x,y) 


^  ~  JVxU^   f       where    f  =  f(7)  only 

The  stream  function  was  defined  to  satisfy  the 
continuity  equation. 


u  =   a^   =/^xU«  d^  ^     =  U^  f,  where   f '  =  di 
ay  d7  ay  d^ 


(2.13) 


au   =   U^   di'  i:?  =  Ucp  f "  y   Ufli 

ax       d7  ax         J  v* 


I 

2x^' 


2x 


7 


V  = 


^^ 


-;^ 


=  "  Ji^xU„   dl  ^7    -  JTu^   (  1  /  2  x   )  f 
ax  d  7  ^x 


=   1/2 


/^ 


(f7   -f) 


a  u   =   Ua, di'  a^   =  U^  f'7  Ua> 


ay 


d-7   dy 


y  X 


(2. 14) 


(2. 15) 


(2. 16) 
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au  =  u^  _u<» 


df"  a?   =     Ua>    f" 


(2. 17) 


^  y^     ^     V  X  d^  ay       ^  x 

Substitute  Eqs .  (2.13)  "  (2.17)  into  Eq .  (2.12),  then 

f£"  +  2  f " '  =  0 

(E.C)    T=0;  £'=0,   f=0 

^   =co  ;  f '  =  1 


The  solution  of  this  Blasius  equation  is  presented  in 
"Boundary  Layer  Theory"  by  Schlichting.   From  the 
transformation  relation. 


JL   = 

X 

(— )  ^ 


C  '  = 


5.0   =   5.0 


JU^  x/y   fR7 
0. 8604 


_/   =   1.7208 

X 


9 


'R. 


=   0.664 


fR, 


0  .  222 /C  U 


CO 


0.664 


^> 


b.   Falkner-Skan  Method 

The  Falkner-Skan  transformation  is  for  2-D, 
axisymmetric  laminar  flow.   The  simplified  Navier-Stokes 
Equation  is: 


u  b  u 

ax 


V  a  u   =   -  1 


ay 


/^ 


IP 
a  X 


^  i^ 


ay" 


(B,C)      y=0;    u=0,   v=0 
y=<»;    u=U(x) 

Take  the  same  t^  as  Blasius'  but  different  with  a  function 
f  +  f(x,'^)  and  follow  the  same  procedure  as  Blasius'  using 
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-  J^  _^   =   U dU 

Z'  dx        dx,   then  the  transformed  equation  is 

f"  +   m  +  1  ff"   +   m 


1  -  (  f '  )""  1  =  X  (  f  _iJ.'   -  f"_iJ_    \ 

•-  -■         ^      ax  ax; 


where  m  is  a  dimensionless  pressure  gradient  parameter 

defined  by: 

m  = X   dU 

U   dx 

(B.C)   ^  =  0;   f  =  0,   f  =  0 

^  =  O);    f  '  =  1 

For  similarity  in  2-D,  laminar  flow,  assume  f  is  a  function 

of.  7)    only  .   Then, 

f"  +   m  +  1  ff"  +   m 


1  -  (f •  ) 


=  0 


(2. 18) 


(B.C)   7=0;   f  =  constant,   f  =  0 
1    =  CO  ;       f  '  =  1 

The  fact  that  m  is  a  constant  leads  to: 

U  =   c  x^ 

where  c  is  also  a  constant. 

In  the  case  of  m  =  0^  ie .  U  is  a  constant, 

Eq .  (2.18)  reduces  to  the  Blasius  equation.   I  f  m  =^  1  which 

means  a  2-D  stagnation  flow,  Eq .  (2.18)  becomes  the  Hiemenz 

equation, 

f"  '  +  ff"  -  (f '  )""  +  1  =  0 

Some  solutions  of  the  Falkner-Skan  equation  for  various 

values  of  m  are  presented  by  Cebeci  and  Bradshaw. 

2 .   Integral  Methods 

a.   Integral  Momentum  Equation 
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For  steady,  2-D  and  incompressible  flow, 
-  1    ^  P   +  .;  a 


u  ^  u   +   V  ^  u   =   - 

ax       ay  1°          d  X 

^  g   +  ^  V  =  0 
^  X     d  y 

From  Eq.  ( 2. 20  )  ,   u  /  a  u  + 

V  a  X 


^  ay" 


(2. 19) 


( 2. 20) 


3  V  N  -  0 


ay 


) 


-  1  IV 


At  the  edge  of  the  boundary  layer,  U  (x)  3  u ( x ) 

a  X 

Then,  Eq .  (2.19)  becomes; 

_a_u  "^  +   a  (uv)   =  U  (  X  )  d  U(x)   +  ))  ±\y_ 
^  X       ay  dx         d':i^ 


Integrate  this  equation  with  respect  to  y,  from  y  =  0  to 

y 


/,  using  7  "^  y/    ^  u 
(  .iu,""  dy   -     U(x)  Jji  dy  -  J  U(x)  dU(x)  dy 


dx 


/^  (2.21) 


Also  we  know,   6     -     \    V-    ~        u   j  dy 

Jz>  ^      U(x)  / 

u(x)   V        u(x)y 


e  = 


dy 


(2.22) 


Substitute  Eq .  (2.22)  into  Eq .  (2.21),  then. 


dU  (x)a   +  /  U(x)   dU(x) 
dx  dx 


r 


w 


r 


This  equation  is  known  as  the  momentum-integral  equation  of 
boundary  layer  theory,  or  as  von  K^rm^in's  integral  equation. 
b.   Pohlhausen's  Method 

Assume  a  polynomial  of  the  fourth  degree  for  the 
velocity  function. 
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u 


U(x) 


J.  3  4- 

f  (A)  =      a  A    +        b  A    +        c  A    +        d   A 


where  A  is  the  dimensionless  distance  from  the  wall, 

^  =   y//   and  0  <  ;i  <  1 

(B.C)   ;\  =  0;    f  =  0, 

^   =    I;         f  =  l,  f'=0,   f"=0 

Also    U(x)  d  U(x)   +  j;  ^  g   =  0   when  ;\  =  0 
dx  a  y-^ 


/   d  U(x) 
y    dx 


-  2b   =  A 


where  J\    is  the  shape  factor. 

then,    f  (A)  =   F(A)   +  il  G(A) 

J     4 
where   F(A  )  =2A-2A  +A 

G(A)  =  1/6  ^  (1  -/^)^ 


Thus  the  boundary  layer  parameters  if  ,    $   and  T^y  can  be 
determined,  if  the  velocity  profile  is  known,  as  follows: 


/*. 


10 


120 


±     =    1   /37 
/      63   I  5 


15 


-'  ) 

144  / 


c.   Thwaites'  Method 

The  integral  momentum  equation  can  be  written 


as : 

dO   +  9  (H  +2)   dU(x)   =  CJ_ 

dx     U(x)  dx        2 

(B.C)   y  =  0;  _£n      =    -  U(x)  K  ,  _£u.   = 
9y^       Q'^  dy 


(2. 23) 


U(x)     L(K) 

e 
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wh( 

sre 

H 

= 

^^ 

C/ 

— 

Tyy 

e 

2 

r 

U"(x 

) 

K 

= 

^ 

dU(x) 

^  L 

e  u(x) 


V        dx 
then   Eq.     (2.23)    can   be    reduced   again, 


U(x)       d_^      =      2      I    L(K)       -      ThCK)    +    2]]    =      F(K)  (2.24) 

y         dx  ^  -^ 

Thwaites  writes  an  expression  for  F(K) 

F(K)   =   0. 45  -  6K 

Then  we  can  write  Eq .  (2.24)  as 

U(x)^    =     0. 45    (  U  (X)   dx 
))  U^(x)  )o 

If  Q    is  calculated  for  a  given  external -velocity 

distribution,  H  and  Cr  can  be  determined  with  the  following 

relations . 

For   0  IK  <_  0.  1 

L  =  0. 22  +  1.57K  -  1.8K 

H  =  2.61  -  3.75K  -  5. 24K 

For    -0.1  <_  K  i  0 

L  =  0.22  +  1.402K  +      0.018K 

0.  107  +  K 

H  =     0.0731   +   2.088 
0. 14  +  K 

3.   Finite-Difference  Methods 

The  finite  difference  methods  are  the  most  flexible, 

practical  and  efficient  methods  to  solve  the  boundary  layer 

equations . 
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The  box  method,  presented  by  Keller  and  Cebeci,  is 
introduced  here  as  the  preferred  finite-difference  method. 
The  momentum  equation  achieved  by  the  Falk.ner-Sk.an 
transformation  can  be  rewritten  in  terms  of  a  first-order 
system  of  PDE ' s .   Then  the  resulting  non-linear  system  is 
linearized  by  Newton's  method.   Finally,  the  block 
elimination  method  is  used  to  solve  the  linearized 
difference  equations  of  the  boundary  layer  problem. 

a.   Box  Method. 

Using  the  new  coordinate  system,  f  and  fi,  the 
transformed  momentum  equation  for  steady,  2-D, 
incompressible  flow  becomes: 

f '  =  u  (2.25a) 

u'  =  V  (2. 25b) 

(bv)'  +  m  +  1   fv  +  m(  1-u  )  =f/u  h  u   -  v  _J_L\ 
2  ^   ^i       ai  / 

(2. 25c) 


n  P- 


Figure  2.1  Net  Rectangle  For  Difference  Approximation 


where   |  =  x  and   t,  =  1  +  ^t/  y 
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A  prime  denotes  differentiation  with  respect  to  "f  . 
The  boundary  conditions  are 

T  =  0;   f ( ^,  0)  =  0,   u(  ^  ,  0)   =0        (2.26) 

We  denote  the  net  points  shown  in  Figure  2.1  by 
.%=0;  f=  f  +Krt    n  =  l,  2 n 


f.  =  °'  7/=  ^/ 


/ 


+  h; 


j  =  1,  2, .  .  .  ,  J 


(2. 27) 


Here  n  and  j  indicate  sequence  numbers. 

With  g-  =   g(  |n,<j)   denoting  the  value  of  any 
quantity  at  the  mesh  point  (  f n,  t  j),  centered  quantities 
can  be  written  as  two-point  averages: 


n-'^. 


r»   n-i 


^       =1/2(^+1  ),  7j.>^      l/2(  7 j  +  7j-l)    (2.28a) 
?^'  =  l/2(  ?y  ^  ?^-).   ^^,=   1/2(J^.^^.,)    (2.28b) 
The  finite  difference  approximations  of  Eqs .  (2.25a) 
and  (2.25b)  for  the  mid  point  of  the  segment  P,  ,  P_^   are: 

(2. 29a) 

( 2. 29b) 

Eq.  (2.25c)  can  be  approximated  similarly  by 
centering  about  the  mid  point  of  the  rectangle   P1,P2,P3,P4 
1-   Centering  Eq .  (2.25c)  about  the  point 

( I  ,  ^) 

Let  the  left-hand  side  of  Eq .  (2.25c)  be  L  and  use 
Eq.  (2.28b) 


n 

— 

f) 
f.f-/ 

= 

n 

n 

^• 

n 

Ui^•-l 

■^^ 

29 


n  n-i 

1/2    (L      +      L         ) 


n-'/i 


L       ^      k,        /  1      k,        /J 


2K, 


(a    )       -    (u    )         -    (fv)         +    f        V      -   V        f    +    (fv) 


By    rearranging. 


(bv ' )  +         m    +    1    ( fv)       +    m 


2 


1      -      (  U      ) 


=    -    L       +  -2 [(a'')     -     (K*) 


-(fv)       +f         v-v         f      +(fv)        J 


(bv')       r     ^,  (fv)       -cx^(u    )       +  o<(v"     f       -    f      V    )       =      R         (2.30) 


where  <X    = 


^.-k. 


R         =    -   L      +  oi     (  fv) 


[< 


n-l 


[(bv)  • 


(u    )       J      -   m" 
+      m    H-    1    (fv)       +    m(l    -    u'^)J, 


fl- ) 


t 


/j'K 


2.       Centering   Eq .     (2.30)    about    the    point 
,     y  .         )       Using    Eq .     ( 2. 29)  , 

r  1  n        n  n        n 

|(bv)'J.=      bj    Vj     -    b».)V^-, 


The    Eq.     (2.30)    becomes 
I  1        rt  n 


h^      (b^   V,-   -  b^-.,   V^,  )    +-<,(fv)^-,;^  -^,(U^)y.;^ 


where 


+  cK    (  f  V  )^'-;£ 


^""^/-x] 


n-i 


-  m 


(2.31) 


^J-'/i  "  l^^i'  ^^^■^?"  "  ^^'  '^^->    ^       ^   nL_L_JL  ^f^V-X"^   m(l-u^)^-.j4.j 


n-i 


The  boundary  conditions  at 


-f- 


are 


f.   =   0- 


u  '   =   0 
o 


V  =   ^' 


(2. 32) 
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b.   Newton's  Method 

For  simplicity,  let  (f-  ,  a-  ^  v^  )   be 
(  £j,  uy,  v^-  )  at  i  ^  ^ 
Then  Eq.  (2.29)  and  2.31)  can  be  written  as 

f^'  -   f/.,   -     M.       (^a  +   ^^'-1   )   =  °  (2.33a) 

^i    ~      ^i'\        ~    -^       (^^'  +   v^^-i   )   =0  (2.33b) 

2 

-t  J. 

h-  (b^  vy  -  b^'-,   v^'.,  )  +   (?<i(fv)y-;i    -  ^^^  ( u   )^'-ji 

+  c^(vj-.^  f/-)i  -  f  J_;i  v^'-y^   )  =  R^-_^       (2.33c) 

n-l 
Here  the  unknowns  are  on  the  left-hand  side  and  R j  -  )i 

involves  only  known  quantities. 

Now,  Newton's  Method  is  applied  to  turn 

Eq.(2.33)  into  a  linear  system  complemented  by  boundary 

conditions  Eq.(2.32)  and  initial  values 


(1  <j  <  J-1) 


The  superscripts  in  parenthesis  represent  the  iteration 

number  as  follows: 

(<;     (o     (h 
t^   r       u^-  ,   v^-  i  =  0,  1,  2,  .  .  . 

where  /f  <<  f,   /u  <<  u,  /v  <<  v. 

Replace  t-i,    m^  ,    v^^  in  Eq .  (2.33)  with  these  expressions. 
Then  Eq. (2.33a)  becomes: 


0 

Co) 

Uo      = 

0 

Vo      = 

v/-' 

<- 

rt-l 

"7  = 

u  . 

^;" 

^7  = 

£J"' 

1 

(O) 

'^r' 
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Thus, 


^^^'  ■   '^^^'-'  "  -^['^^Z  ^  '^^^■-'  ]  =    ^'^i 


(2. 34a) 


(1^ 
where    ( i^  i  )  «-    =    £  ^■- 1 


(I) 


(I) 


^^'     ^      h^'  ^r^x 


From   Eq.     (2. 33b) , 


(2. 34b) 


(i;  (I)  (i; 

where    (rj)^-,,=      u,>.,    -   u^-      +    hy   v^..^,^ 


From   Eq.     (2. 33c) 


-<  ''^  /     ii;  (t)         J.      It)  f  (n  ^       ^     ID 

h    ■    (     h^^v^    -    b^-.,     ^  V^L,      )     +      (?<,  <7  (  fv)y.  ^^    -     o<j^(u     )j^>i 


n-i      ^      <i; 


here     /(£v,;:^=l/2     [f^fVvV'     .    .f  f  ^^  ^      i-.,i^^,     .      «;-/£;(',] 


ff      (    U  )4"     |y        =  U     ■  </      U    V  +  U  V  /7     U    V 


/u)  (I)  f.      in 


^  ^/-)i 


1/2  ( /v^-u/v;:;    ) 


then    (s,)j,'  /v^- +    (Sx)^Vv^'.,     +    ( Sj  )h"  /  £."    +    (s^)-'/fv.,    +    (Sy)//u/ 
+    (s.  ).^     /u/  .  =       (r,).^  ^  (2.34c) 


(r.), 


where    (r^)^^    =   R^-.^^     -      [h^-Cb^-   v^-    "   by.,      v^.,       +    -^.(£v)/.y^ 


<^         f 


f\-\ 


d-'/i^ 


^  ^         ~J~        ^         "2~        ^■^- 
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(s^),- 


2  2 


-  o/a  U 


1'/ 


the  boundary  conditions  are: 
/f(7   =0     /uo  =   0 


i-uj  = 


0 


(2. 35) 


c.   Block  Elimination  Method 

Eqs.(2.34)  and  (2.35)  are  the  linearized 
difference  equations  of  the  momentum  equation  for  external 
flows  which  has  a  block,  tridiagonal  structure.   We  can  solve 
these  equations  by  means  of  the  Block  Elimination  Method  as 
discussed  by  Keller  (1974). 

Let  us  define  the  three-dimensional  vectors,  /v 
and  r^,  to  express  the  system  in  matrix-vector  form, 

0 


/,  = 


0  <  j  <  J 


^o    = 


^^5     K 


'^ 


(r, 

'^ 

(r;, 

'J 

"^J 

'J 

1  <  j  <  J-1  ( 2. 36) 
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(r,  )j 

0 


J  =  J 


and  the  3x3  matrices,   Aj  ,  Bj-  ,  C.^ 


A.   = 


10  0 
0  10 
0    -1 


-h,  /2 


1         -h^W2         0 
(Sj  )^-     (S^-y        (S,  )y 
0  -1  -h^'^,/2 


1    <    j     <    J-1 


Aj   = 


Bi    = 


1  -hj/2  0 

0  1  0 

1  -h//2  0 
0              0  0 


j  =  J 


1  <  j   <  J 


(2. 37) 


0 

0 

0 

Cj  = 

0 
0 

0 
1 

0 

0  <.  j  <.  J-1 

Then  Ec 

IS. 

(2 

34) 

and  (2.35 

) 

can  be  written  as 

/^  /=     R 


(2.38) 
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where 


P   = 


B,  A,  C, 


^i    A  J  q 


Bj..  Aj-,  Cj., 
B7    Aj 


/  = 


f^  = 


i^. 


Let  us  factorize  the  matrix  ^    to  solve  Eq .  (2.38) 

/?  =   xy  (2.  39) 


where 

I 


P> 


X,  I 


x^-  I 


X 


y  c 


y,     c, 


^3-  ""i 


Vj-.  ^J-. 


Here  I  is  the  3x3  identity  matrix. 

X  ^-   and  y-   are  also  3x3  matrices. 
According  to  Eq.(2.39),  we  can  find 

yo  =   Ao 

X3  y^-.,  =   B^"   (  j  =  1,  2,  .  .  .  ,  J) 

^i    =  ^;  -  ^/  ^/-.    ^^  =  ^'2..  .,  J) 


(2. 40a) 
( 2. 40b) 
(2.40c) 


and  the  matrix  xv  has  the  same  structure  as  that  of  the 
matrix  B*',  Therefore,  x^-  has  the  elements  like  this. 


^J  = 


(x  „  )y    (X  ,i  )^'    (X,j  )^' 

(Xai  )-•   (X;u  )v   (x^j  ).- 

0   '     0   ^     0   ^ 
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and  y • can  be  denoted  as: 


(yn 
(yai 

0 


From  Eq.  (2. 40a) , 

(y  K  )o  =1    (y.i 

From  Eq.  (2. 40b) 

(  X  M  )  I  =  -1    ix  ,± 
(X  ^,  ),   =  (S4),  (X  ij^ 

From  Eq.  (2. 40c) 


3   (y.j.  U     ^Y,3   )y 

-1      -h,4,  /2^ 


=  0     (y,3  )o  =  0 
=  1     (y  . )   =  0 

"'  a  J  o 


=   -1/2  h,  (x,3  ),  =  0 

=  (s^),  +  (Xi;,),   (x^3),=  -2   (s^),/^, 


"^    -         '^     2  2 


Then  we  can  compute  the  elements  of  x   with  Eqs.  (2.37)  and 
(  2.  40b)  for  1  <_  j  <_  J, 
(x„  ).'  = 


(  X  ,^  )^^  =  - 


Ya  ^    (_  Z 

L_|(^)^  -    (x„  ),^    [(y,.  v.,  ^     -    (y,3V-.]  1 

(Xi,  )^^   =  1_  j(Sa)/(yai  ^^'-^     +  (  S4  ),"  (  y  ;,3  )^>- ,  + 

Yo   (  A*  r(s^V(y^  )^i,-  (s^)^>(y^,  )j>_,  j 

(J2i(s,)^^  -  (s^)^v  +  (xa,  )^^  [(y,^  )^i,-_h^"  (y/J^^,]  I 


-^v=^ 


(x^j)^-  =   (x^,  )^^  (y,^  )3'.,   +  (x^)y(y^  )j._,   -  (s^  )^^ 
^0=  ^y/3^3-1  ^y-i  V-.   -  (^:^3)j^.,  (y.,  )^\y   -  -, 
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y,  =  ^  (y^V.,  -  <y.,',- 

Let  y  /  =  z,  then  xz  =  R 


Thus,  z 


o    =  ^o 


(2.41) 
(2.42) 


where 


^^"  ^r 


(1  <  j  <  J) 


(z,  )^. 

(  Z,  )  ' 


(0  <  j  <  J) 


(23)0  =  (^^5)0 


From  Eq.  (2.42)    for  j  =  0 

(z,  )c  =  (^.  )o  (2,)^  =  (r,  )^ 
and  for  1  <_  J  <_  J 

(z,)^-  =  (r,)^-  -  (x„  )^«(z,V-,  -  (x,^  )^>  (z^)^^,  -  (x,j  )y(Zj)^'_, 

(Z;^)y     =    (ri)y     -     (X;,,   )^>   (Z,   )^l,     -     (X^   )^^   (Z^  )^:,      -    (X^^   )^'  (Z3  )^'., 

(Zj)/  =  (rj)^^ 
From  Eq.  (2.41), 

Yj  /j  =  z-j-  (2.43a) 

y^V/  =  ^i  -   ^^V^Vi    (0  1  3  i  J-1)       (2.43b) 
Then  the  vectors    can  be  calculated  with  Eq .  (2.43).   The 
three  components  of   ,  for  j  =  0, 1, 2, . . . J-1,  are: 


/U, 


2  ' 


/v^-  =  I  ((y.,),«[(z,)^^  +  ^'^y- V]  -  (y-v^^'V  -  ^'^y.,V^^uV  ] 

/w-    =^_^      [(z,)-     -    (y^^   )^^     /-u^^    -    (y,^)^-   /v^>] 
(Yn  V  -^ 

where    e,  =    (z^)^-    -    /  u^^^,    +    _h£4,     /"v^-^^ 
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-   (y^,  )/  (y^j  )/  +  (ya3  )/  (y  n  V 

and  for  j  =  J 

/vj  =     e^  (y;,,  )j   -  ej    (y„  )cr 


/w,  = 


^y/3  ^cr^y-.  V     -  (yx>)j(y  ..  v 

where   e^  =  (z,)j  -  (y,^  )j-/u3- 
These  calculations  are  stopped  when 

l/v/''  I    <    e 

where   v(o)  is  the  wall  shear  parameter 
^   is  a  prescribed  value. 


F.   TURBULENT  BOUNDARY  LAYER  CALCULATIONS 

Turbulent  fluid  motion  is  an  irregular  condition  of  flow 
in  which  the  various  quantities  show  a  random  variation  with 
time  and  space.   Therefore,  turbulence  is  characterized  by 
random  and  chaotic  motion  of  fluid  particles. 

The  velocity  varies  randomly  at  any  point  in  a  turbulent 
fluid.   The  velocity  components  in  a  three-dimensional 
turbulent  fluid  are: 

u   =   u   +   u ' 

V   =   V   +   V  ' 

w   =   w   +   w ' 
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where  u,  v,  w  represent  the  mean  velocities 
u   =      1 \    u   dt,  and  so  on. 


t.  -  t,   . 


u',  v',  w*  represent  the  amount  of  fluctuation  of  the 
instantaneous  velocities  from  the  mean  velocities. 


Obviously,  we  can  see  that 
1      I  u'   dt   =   u' 


o,    etc 


'     i. 


Similarly,    the    following    quantities    can   be    defined, 
1  I  u'"^  dt      =      u'""    ,    etc. 


tjL    -      t,   -^^ 


1  u'v'      dt      =   u'v',      etc. 

Even  though  the  flow  is  turbulent,  the  time-dependent 
momentum  equations  are  valid.   However,  it  is  impossible  to 
solve  the  momentum  equations  for  turbulent  flow  because  the 
fluctuations  are  random  and  chaotic.   Thus,  Reynolds 
modified  the  momentum  equations  by  introducing  the  mean 

a 

value  and  fluctating  values  of  the  flow  quantities  and  by 
assuming  the  fluctuations  to  be  continuous  functions  of  time 
and  space.   From  the  momentum  equations. 


<»u    +    u   du    +    V    ^u    +    w   <^u 


-1      ^P      +>>/_iji      +_/u      +_^  \ 
1^      ^x  i  c?x^         dy*         <^zV 


^t  ax  ay  a- z  /3      <^x  V  <?x-         dy 

By   multiplying    the   continuity   equations    with    u 

=   0 


u  /_^_u      +     ^  V       +      ^  w  A 
\  d  X  ay  ^  z   / 
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Sum    these    two    equations    and    rearrange,    then, 
3u      =    -    1    ^P      +      a 


^  t 


C  ax        ,)x 


V<ju    -    u^V      ^    I  l^ifa       -    uv")    +      ^  /  ^^\i 


ay     ay 


<^z  I    <)  z 


UW     j 


Here    take    the    mean    value    of    each    term   using    u'    =    0, 


^(u    -I-    u '  )         =_i_u      +_iu.'       =     <?u 

at  at        at  at 


-  _!_     a(P  -t-  P' )     =     ^    ap 
/**  ax  /^     a  X 


^(  u    4-    u '  )         -    (u    +    u')"^       =^au      -    (u-^    +    2    uu'    +    u) 


ax 


ax 
=  PaQ 


v> 


^ 


a  (u    +    u '  )         =  V'^u      -      uv      -   u'v' 


a  y 


ay 


acu  +  u' )       =  ^aii 


UW       -    u ' w 


(2.45) 


^  z  a  z 

Substitute    Eq .     (2.45)     into    Eq .     (2.44)    and    rearrange.       Then, 


/O  ( "iiL    "*"     ^  au     +     V  a_5     +     w  au    \ 
I    \  dt  ax  ay  az     / 


/" 


=  -_AP_+^  c^u"  /^(  a  u'  +     au' v'     +  au'w'   j 

ax  ax  ^  y  dzM2.  46a) 

By    similar    procedures    for    the    other    directions, 
(av      +uav^      +v    c>\/      +    way      \ 


at 


^A 


^y 


dZ 


=    -_£_P 

ay 


+/* ^v  -^/  a V 'u'   +  j|_v'    +    a  V ' w '  ^ 
'^l     a  X         ay  d  z    i 


(2.46b) 


r 


(  a  w     +  Q     a  w     +  V      a  w     +  w      a  w 


at  d  X 

=  -    ap 


a  y 


a  z 


) 


az 


+/r^*w  -x7/_iw_|_u'   +    ^  w ' v'     +      a  W  '    \ 
'(ax  ay  az^ 


(2.46c) 


and  from  the  continuity  equation, 

3  u   +    a  V   +   a  w   =  0 
ax     d  y     a  z 


(2.47) 
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In  laminar  flow,  we  have  4  unknowns  and  4  equations. 
But  in  turbulent  flow,  we  have  10  unknowns  and  only  4 
equations  as  shown  at  Eqs.(2.46),  (2.47).   This  is  the 
reason  why  we  need  the  turbulence  models  which  were 
discussed  in  Section  C. 


G.   SEPARATION 

In  some  cases,  the  boundary  layer  thickness  increases 
considerably  in  the  downstream  direction  and  the  flow  in  the 
boundary  layer  reverses  its  direction.   The  change  in 
direction  causes  the  decelerated  fluid  particles  to  move 
outwards,  which  means  that  the  boundary  layer  is  separated 
from  the  wall.   This  phenomenon,  boundary  layer  separation, 
is  always  related  to  the  formation  of  vortices  and  to  large 
energy  losses  in  the  wake  of  the  body. 

Let  us  consider  the  simplified  boundary  layer  equation 
in  order  to  investigate  two-dimensional  separation; 
u   3  u   +  V  ^  u   =  -  1    ^P   +  \ld_a 

<^x     ay     /o       ax    ay^ 

and  define  the  separation  as 
/  ^  u  ^     =  0 

Since   u=v=Oaty=0, 
^)  /  '^^I  \     =1     ^P 
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^^^ly-^o  ^  ^^ 


The  velocity  profile  in  the  boundary  layer  always  has  an 
inflection  point  in  the  region  of  decelerated  flow  and  the 
velocity  profile  at  separation  point  must  have  an  inflection 
point.   Thus,  separation  can  occur  only  when  the  flow  is 
decelerated.   In  two-dimensional  separation,  a  bubble  of 
fluid  which  has  low  velocity  is  always  formed.   This  bubble 
is  often  unsteady  and  distinguished  by  interior  streamlines 
which  are  closed  loops  or  which  extend  to  infinity. 
There  are  two  severe  problems  in  boundary  layer 
calculations  when  separation  occurs. 

1.  The  Goldstein  singularity  at  the  separation  point  in 
direct  boundary  layer  calculations. 

2.  Numerical  problems  downstream  of  the  separation  point. 
If  a  boundary  condition  prescribes  the  pressure 

gradient,  then  the  boundary  layer  methods  suffer  from  a 
singularity  due  to  separation.   Moreover,  the  singularity 
causes  the  numerical  breakdown  in  direct  boundary  layer 
methods  near  the  separation  point.   This  Goldstein 
singularity  can  be  overcome  by  using  the  displacement 
thickness  or  the  wall  shear  stress,  instead  of  the  pressure 
gradient,  for  the  boundary  condition,  which  makes  it 
possible  to  integrate  the  boundary  layer  equations  through 
the  separation  point.   Also,  the  full  Navier-Stokes 
equations  exhibit  no  singular  behavior. 
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On  the  other  hand,  the  boundary  layer  equations  lead  to 
a  numerical  instability  in  the  region  o£  reversed  flow.   The 
FLARE  approximation  is  the  most  common  method  to  overcome 
this  instability.   The  momentum  transport  term 
u  ^u/<?x   is  deleted  where  u  is  smaller  than  zero.   But  the 
accuracy  of  this  approximation  decreases  as  the  region  of 
reversed  flow  increases.   Therefore,  it  is  necessary  to 
introduce  an  upstream  influence  in  that  region.   This  is 
incorporated  by  the  downstream  upstream  iteration  procedure 
(DUIT)  which  consists  of  a  sequence  of  alternating  up-  and 
downstream  sweeps  with  the  momentum  transport  term   u  au/c?x. 
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III.   INTERACTION  METHOD 


A.   INTRODUCTION 

If  the  flow  remains  attached  and  the  Reynolds  number  is 
high,  the  pressure  distribution  and  the  overall  lift  force 
can  be  obtained  from  a  potential  flow  solution. 
Conventional  boundary  layer  methods  provide  additional 
information  about  the  skin  friction  distribution  and  the 
overall  drag  force.   However,  if  the  flow  separates,  no 
information  is  available  for  regions  downstream  of  the 
separation  point.   Therefore,  the  overall  forces  cannot  be 
obtained . 

For  this  reason,  interaction  methods  are  introduced  to 
overcome  this  problem.   They  provide  a  special  coupling 
between  the  inner  viscous  and  the  outer  inviscid  flows,  and 
can  be  classified  as  follows: 

1.  Weak  interaction  methods 

a.  Direct  method 

b.  Inverse  method 

c.  Semi-inverse  method 

2.  Strong  interaction  method 

The  weak  interactions  provide  only  a  loose  coupling 
between  viscous  and  inviscid  regions,  i.e.  two  different 
regions  are  treated  alternately.   The  viscous  flow  solver 
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deals  with  the  flow  in  the  viscous  region  and  yields  the 
boundary  condition  of  the  inviscid  region,  while  the 
inviscid  flow  solver  deals  with  the  flow  in  the  inviscid 
region  and  yields  the  boundary  condition  of  the  viscous 
region.   The  exchange  of  information  through  the  boundary 
condition  is  slow,  but  regions  such  as  those  near  separation 
and  trailing  edges  require  fast  and  direct  coupling  between 
viscous  and  inviscid  region.   While  the  strong  interaction 
method  treats  displacement  thickness  and  external  velocity 
simultaneously,  the  weak  interaction  methods  process  one  of 
these  quantities  as  input  and  the  other  as  output. 

1.  Inviscid  flow  methods 

a.  Direct  boundary  conditions 

(1)  Prescription  of  the  airfoil  shape 

(2)  Zero  normal  velocity  at  the  surface 

b.  Inverse  boundary  conditions 
Prescription  of  a  velocity  distribution  for 
the  unknown  airfoil  shape 

2.  Viscous  flow  methods 

a.   Direct  boundary  conditions 

(1)   No  slip  condition  requiring  zero  normal 
and  zero  tangential  velocity  at  the 
surface . 
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(2)   Prescription  of  the  external  velocity 
i.e.  a-component  of  velocity  at  the 
edge  o£  boundary  layer 

b.  Inverse  boundary  conditions 

(1)  No  slip  condition 

(2)  Prescription  of  the  displacement 
thickness 

c.  Boundary  conditions  for  simultaneous 
interaction 

(1)  No  slip  condition 

(2)  Prescription  of  a  linear  combination  of 
displacement  thickness  and  external 
velocity 

B.   WEAK  INTERACTION  METHODS 
1 .   Direct  Method 

A  direct  inviscid  flow  solver  is  combined  with  a 
direct  viscous  flow  solver  (see  Figure  3.2).   The  boundary 
conditions  are: 

u  (x,o)   =  0    v(x,o)    =0  (3.2) 

u  (x,  y  )   =   Ue  ( X) 
This  method  is  terminated  at  the  point  of  vanishing 
skin  friction.   In  the  direct  scheme,  the  pressure  is 
calculated  from  the  inviscid  region,  but  the  displacement 
thickness  is  determined  from  the  viscous  region.   Because  of 
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this  phenomenon  and  Goldstein  singularity,  this  method  is 
not  appropriate  for  flows  with  strong  interference  effects 
between  viscous  and  inviscid  regions. 

2 .  Inverse  Method 

This  method  consists  of  an  inverse  inviscid  and  an 

inverse  viscous  flow  solver  (see  Figure  3.3).   The  boundary 

conditions  of  the  inverse  boundary  method  are: 

u  (x,o)   =  0  ,      V(x,o)   =0  (3.2) 

//»  (x,y^  )  =   u(x,  yg)   [Ye  -  J'*(x)3 

where  ^  is  the  stream  function. 

The  roles  of  displacement  thickness  and  external 

velocity  distribution  are  exchanged  in  this  method.   The 

troubles  related  to  the  Goldstein  singularity  can  be 

overcome,  but  the  whole  procedure  takes  a  long  time  due  to 

very  slow  convergence.   Thus,  this  method  can  be  applied  to 

the  regions  of  separated  flow  only  and  needs  severe 

under -relaxation. 

3 .  Semi-inverse  Method 

A  direct  inviscid  flow  solver  is  combined  with  an 
inverse  viscous  flow  solver  (see  Figure  3.4).   The  input  is 
the  displacement  thickness  and  the  output  is  the  external 
velocity  distribution  for  both  solvers.   The  two  external 
velocity  distributions  are  combined  and  then  an  updated 
displacement  thickness  is  obtained  through  a  relaxation 
procedure.   A  formula  for  satisfactory  convergence  is: 
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^  new    (x)       =     <7 


old    (X) 


1      +      w   / 


UeKAJ 


(3.  3) 


where  w  is  the  relaxation  parameter. 

The  numerical  weaknesses  can  be  overcome,  but  the 
coupling  is  still  loose. 


C.   STRONG  INTERACTION  METHOD 

The  simultaneous  method  solves  the  boundary  layer 
equations  subject  to  an  interaction  law  (see  Figure  3.5). 
Viscous  displacement  effects  are  allowed  to  cause 
substantial  changes  in  the  external  velocity  distribution. 
Since  both  displacement  thickness  and  external  velocity  are 
treated  as  unknowns,  one  more  additional  relation,  the  so 
called  interaction  law,  is  needed.   Thus  the  boundary 
conditions  are: 

u(x,o)   =   0      v(x,o)    =   0 

/A(x,y^)    =   u(x,yg)   [y^-   /(x)]  (3.4) 

u  (x,yg)    =   u/(x)   +  l/ir  [ d__[u(|,y<f)  /  ii)]       di 


X  -I 


The  last  equation  represents  the  interaction  law. 

Here,  the  external  velocity  for  the  boundary  conditions 
can  be  written  as: 

U  (x,ye)   =  Ue   ^x)   =   U^  (x)    +  ^  Ue  (x)       (3.6) 

o 

where   U^  (x)  is  due  to  inviscid  flow  past  the  airfoil 

^Ug     (x)  is  the  perturbation  velocity  due  to  the 
displacement  effect  of  a  boundary  layer. 
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To  obtain  ^U©.  ( x ) ,  the  blowing  velocity  concept  is  used.   Let 
us  consider  an  airfoil  along  which  sources  are  distributed 
(see  Figure  3.1).   This  surface  distribution  of  sources 
represents  the  effect  of  boundary  layers.   There  are  two 
ways  to  predict  the  displacement  effect  of  a  boundary  layer 
on  the  outer  inviscid  flow: 

1.  Compute  the  inviscid  flow  past  a  displacement 
body 

2.  Replace  the  condition  of  zero  normal  velocity  on  the 
surface  with  a  condition  of  prescribed  blowing 
velocity  at  the  surface 


Y   T 


> 


Jk  U^(x) 

I   V  ( X,  o )   =     2 


^ ^      X 

Figure   3.1   Blowing  Velocity  Concept 


The  streamlines  are  displaced  away  from  the  surface  by 
the  distributed  sources  which  eject  the  fluid  at  the 
surface.   Then  the  virtual  displacement  body  becomes  a 
streamline  and  the  flow  tangency  condition  is: 


Utf(x)        Jx 


(3.6) 


From  the  thin  airfoil  approximation,  the  displacement 
thickness  is  assumed  to  be  so  small  that  u-component  of 
velocity  do  not  change  across  the  layer  and  the  airfoil  in 
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this  connection  can  be  represented  by  a  straight  line. 
Therefore,  the  blowing  velocity  v(x,o)  is  equal  to  half  of 
the  source  strength  <r(x). 


V  (x,o)   =  1/2  (T  (X)  (3.7) 

/^         r 

v(x,o)   =   v(x,<7  )    -    ^v    dy 

Jo  <^y 

=   Ug  d/_    +   /^dLM 
dx  dx 

=      6      iVe  /    )  (3.8) 

dx 

Thus,    d  iUeS   )    =   1/2  (T  ( x )  (3.9) 

dx 

Now,  the  perturbation  velocity  due  to  the  displacement 

effect  can  be  written  as: 

4Ue(x)   = 1     (  '■tfjjj.   d^  (3.10) 


2-r   ;^  X  -J 

Finally,  Eq .  (3.6)  becomes: 

Ug(x)   =   Ue  (x)   + ^  (   (r(>)   d^  (3.11) 

27r   J^^x  -> 

Eq .  (3.11)  represents  the  interaction  law  in  usable  form  and 

the  integral  on  the  right  hand  side  is  known  as  the  Hilbert 

integral . 

In  the  simultaneous  method,  the  inviscid  flow  solver 

provides  an  initial  external  velocity  distribution,  but  the 

inviscid  flow  solver  is  not  incorporated  in  the  overall 

iteration  process.   Thus,  the  viscous  flow  solver  needs  the 

interaction  law  to  obtain  the  necessary  information  about 

the  inviscid  region. 
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VISCOUS  FLOW  SOLVER 


(Direct  Boundary  Conditions) 
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dis tr  ibut ion 
Output:  Displacement  thickness 
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l</(x) 


No 


Figure   3.2    Direct  Method 
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VISCOUS  FLOW  SOLVER 


(Inverse  Boundary  Conditions) 
Input:   Displacement  thickness 

distr  ibut ion 
Output:  External  velocity 

distribution 


4 


U^  (X) 


Air  fo  il 
Geometry 


INVISCID  FLOW  SOLVER 


(Inverse  Boundary  Conditions) 
Input:   External  velocity 

distribution 
Output:  Displacement  body 

shape 


No 


Figure   3.3    Inverse  Method 
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INVISCID  FLOW  SOLVER 


(Direct  Boundary 
Conditions ) 
Input:  Updated  shape  of 

displacement  body 
Output : External  velocity 
dis tr  ibution 


VISCOUS  FLOW  SOLVER 

(Inverse  Boundary 
Condition) 
Input:  Displacement  thick-  |^ 

ness  distribution 
Output:  External  velocity 
distribution 


RELAXATION 


Update  displacement  thickness 
distribution  thickness  distribution 
and  the  two  external  velocity 
distributions . 


Figure  3.4   Semi-Inverse  Method 
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(Direct  Boundary  Conditions) 
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UgCx) 


VISCOUS  FLOW  SOLVER 
and  INTERACTION  LAW 


(Displacement  thickness  and  external 
velocity  are  treated  as  unknowns) 
Input:   Data  obtained  from  previous 
sweep  and  initial  inviscid 
solution 
Output:   Displacement  thickness 

distribution  and  external 
velocity  distribution 


d   (X) 


4 


OVER 
RELAXATION 

(X) 


OF 


No 


Figure   3.5    Simultaneous  Method 
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IV.   DESCRIPTION  OF  CEBECI ' S  INTERACTIVE  PROGRAM 

A.   INPUT  DESCRIPTION 
1 ,   Introduction 

The  Reynolds  number  is  the  most  powerful  parameter 
which  can  affect  the  flow.   Wholly  different  flows  can 
result  from  different  Reynolds  numbers  (e.g.  6.0E  +  05 
against  0.28E  +  06).   Even  though  a  flow  with  high  Reynolds 
number  behaves  well  and  does  not  separate,  it  may  have  an 
extensive  flow  separation  at  a  low  Reynolds  number. 
Numerical  breakdown  in  boundary  layer  computations  is  partly 
due  to  this  extensive  flow  separation  at  low  Reynolds 
numbers.   To  avoid  unreal istical ly  large  regions  of  flow 
separation,  it  is  necessary  to  make  some  changes  in  the 
turbulence  model  for  transitional  flow.   Reduction  in 
transitional  region  to  a  proper  level  is  a  method  to 
decrease  numerical  problems  in  computation.   By  doing  this, 
more  reasonable  results  can  be  obtained  in  low  Reynolds 
number  flows.   Also,  the  computation  at  low  Reynolds  numbers 
needs  more  iterations  to  converge  than  at  high  Reynolds 
number,  and  the  lift  coefficient  increases  with  Reynolds 
number  at  the  same  angle  of  attack  because  low  Reynolds 
number  flows  exhibit  a  stronger  displacement  effect  than 
high  Reynolds  number  flows. 
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The  results  may  be  affected  by  transition  location, 
The  transition  location  can  be  either  fixed  by  the  user  or 
computed  from  the  following  empirical  formula  given  in 

Cebeci  and  Bradshaw  (1977). 

0.4  6 

Rq  =   1.174  I   1  +   22400  ^  R^  (4.1) 


,74  (  I    +       22400  ^ 


X 

Laminar  flow  separation,  however,  may  occur  upstream 

of  the  transition  location  predicted  by  the  above  equation. 

In  this  case,  it  is  assumed  that  the  onset  of  transition 

corresponds  to  the  point  of  laminar  separation,  and  the 

following  message  will  be  issued  in  the  output; 

TRANSITION  LOCATION  HAS  BEEN  RECOMPUTED  AT  THE 
POINT  OF  LAMINAR  SEPARATION 

The  results  of  this  program  agree  well  with  the 
experimental  results  up  to  stall  at  high  Reynolds  number 
flows.   However,  low  Reynolds  number  flows  may  not  agree 
well  or  experience  numerical  breakdown  close  to  stall  with 
the  following  messages: 

at  the  very  top  of  output, 
+  IFY  002  I   STOP  1 

or  at  the  very  bottom  of  output, 

IFY  207  I  VFNTH:   PROGRAM  INTERRUPT  (Z)- 
FLOATING-POINT  EXCEPTION  OVERFLOW 

IFY  259   I  STNCT:  /ARG/  =  /argument/, 

(HEX  =  hexadecimal),  APPROACHES  SINGULARITY 

or  MAXIMUM  NUMBER  OF  ITERATIONS  EXCEEDED 

Thus,  the  range  of  angle  of  attack  should  be  considered 

carefully. 
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If  the  user  requires  many  sweeps  with  IPRNT  =  0  or 
2,  then  the  user  needs  to  choose  class  "G"  or  "J"  for  enough 
running  time  and  add  additional  command  ",  LINES  =  (m)"  just 
after  nnnnP  on  the  following  line  to  make  enough  space  for 
pr  int : 

11*  MAIN  ORG  =  NPGVMl   nnnnP     ,  LINES  =  (m) 

where  nnnn  is  user's  ID  number. 

m  is  the  number  of  output  lines  in  thousand 

(Generally  10 js  enough) 

If  the  user  does  not  do  that,  the  following  error 

message  will  be  issued  at  the  very  top  of  the  output  and 

printing  will  be  stopped  abnormally: 

lAT  1600  JOB  number  (user's  job  number)  LINES  EXCEEDED 
lEF  4501   user's  job  number  -  ABEND  S722  UOOOO 

2 .   Detailed  Input  Data  Description 

The  input  to  the  computer  program  consists  of  1 

title  line,  3  control  lines  and  airfoil  coordinate.   User 

must  follow  the  data  format  for  specified  column  and  type. 


I  Title  Line   I  FORMAT  (18  A4) 


This  line  provides  any  description  as  desired  with 
any  acceptable  machine  characters. 


Control  Line  1  I       FORMAT  (6  15) 
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Column 

Name 

5 

ITF(l) 

10 

ITF(2) 

15 


IRSTRT 


Explanation 

Transition  flag  for  the 

lower  surface 

Transition  flag  for  the 

upper  surface 

=  1   Point  of  transition 
has  to  be  specified 
by  user.   This  option 
shall  be  used  if  ex- 
perimental results 
are  available  or  to 
avoid  oscillations  of 
the  computed  transition 
points.   In  the  latter 
case,  fix  transition  at 
the  most  upstream 
occurrence  of  computed 
transition . 

=  4   Point  of  transition 
will  be  calculated 
according  to  Eq .  (4.1). 
No  input  is  necessary 
for  XTRL  and  XTRU.   The 
computed  point  of 
transition  will  be 
redefined,  if  laminar 
separation  takes  place 
upstream  of  the 
location  predicted  by 
Eq.  (4.1) 

Boundary  layer  restart 
flag. 


Read  Starting  Solution 


Store  Final 
solution  or  Unit  2 


=  0 
=  1 


No 
Yes 


Yes 
Yes 


Leave  0  for  general  use 


20 


IGLMAX 


Number  of  sweeps. 


Low  Reynolds  number 
flows  need  more  sweeps 
to  converge. 


25 


IPRNT 


Print  flag  to  control 
output  print. 
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=  -1    Summary  print  for  the 
last  two  sweeps 

=  0     Summary  print  for  each 
sweep . 

=  2     Whole  print  for  all 
sweeps . 


Control  Line  2 
Column 


FORMAT  (4  ElO.O) 


1-10 
11  -20 

21  -30 


Name 

RL 
XTRL 

XTRU 


31  -40 


ALPO 


Explanation 

Chord  Reynolds  number 

Fixed  transition  location 
for  lower  surface 

Fixed  transition  location 
for  upper  surface  when  the 
stagnation  point  is  above 
the  leading  edge,  XTRL  may 
be  positive.   When  the 
stagnation  point  is  beneath 
the  leading  edge,  XTRU  may 
be  negative. 

Angle  of  attack  in  degree 


Control  Line  3 
Column 


1  -5 


Name 


MPTS 


FORMAT  (15) 

Explanation 

Number  of  airfoil 
coordinates . 


Coordinate  Data  Line 


FORMAT  (2F  10.0) 


Input  dimensionless  airfoil  coordinates  as  two  columns  in 
one  row  (  x/c,  y/c) .   The  order  is  trailing  edge  — slower 

surface  — >   leading  edge  >   upper  surface  ►  trailing 

edge.   Thus,  trailing  edge  is  input  twice. 
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B.       OUTPUT    DESCRIPTION 

1 .  Introduction 

To  take  the  appropriate  results  from  the  output, 
user  must  check,  the  convergence.   This  can  be  done  by 
comparing  the  convergence  indicators,  lift  coefficient  and 
displacement  thickness  at  trailing  edge,  when  the 
computation  is  completed  successfully  for  the  given  sweeps. 
If  the  convergence  indicators  show  steady  values  over  some 
sweeps  i.e.  each  difference  of  the  convergence  indicators  is 
less  than  1%,  the  results  are  considered  as  converged. 

Sometimes,  the  user  may  experience  failure  in 
computations  at  low  Reynolds  numbers  and  high  angles  of 
attack  due  to  numerical  breakdowns.   These  breakdowns  take 
place  when  Newton  iteration  does  not  converge  near  to  the 
trailing  edge  on  the  upper  surface  or  the  computation  is 
terminated  by  Fortran  error  with  this  message: 

IFY  2511  SSQRT;   ARG  =  argument,  LESS  THAN  ZERO 

Also,  extensive  flow  separation  cause  numerical 
breakdowns  of  the  boundary  layer  computation.   These 
unreal istically  large  regions  of  separation  at  low  Reynolds 
numbers  can  be  reduced  by  changing  the  transitional  flow 
model  . 

2 .  Detailed  Output  Data  Description 

The  output  of  the  computer  program  for  IPRNT  =  -1 
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and  0  can  be  divided  into  three  parts  as  follows: 

1.  Input  data  and  inviscid  lift  coefficient 

2.  Alternating  boundary  layer  parameters  for  each 
surface  and  wake. 

3.  Inviscid  and  viscous  pressure  distributions 

The  following  is  the  list  of  computed  boundary  layer 

parameters  in  the  order  in  which  they  appear  in  the  output 

pr  int . 

Name         Name 

in  print     in  program     Definition     Explanation 

X(NX)        X(I)  1_        Dimensionless  sur- 

c         face  distance  from 

the  stagnation  point 

X/C(NX)      XC(I)  X_        Dimensionless  chord- 

c         wise  distance  from 
leading  edge 

V(1,NX)      VWL(I)        /  d_l_f  \      Dimensionless  wall 

^  ^  '*i^   /w  shear  stress  para- 

meter . 

CF  CF(I)  2  Tw  Local  skin  friction 

/OUg^  coefficient 

DELST       DLS(I)  </  Dimensionless  dis- 

c  placement  thickness. 

dy 


THETA        THETA(I)         _6_        Dimensionless 

c         momentum  thickness 


'00 


0  = 


Q 


u 


u_/'l y—")  dy 

e    V        Ue     I 


UE(NX)       UE(I)  Uel        Dimensionless  exter- 

Ua>        nal  velocity 

computed  from  the 
inviscid  method  with 
viscous  effect. 
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W(NP,NX) 


WNP(I) 


UeV 


D(NX) 


DB(NX) 


IT 


NP 


D(KX) 


DB(KX) 


Uev 


i^/RL 


>co 


Dimenslonless  exter- 
nal velocity 
computed  from  the 
interactive  boundary 
layer  method. 

Product  of  velocity 
and  displacement 
thickness  in  dimen- 
slonless form  for 
the  current  sweep 
( D ) ,  from  the  pre- 
vious sweep  ( DB ) . 
D  and  DB  are  printed 
in  the  output  as  up- 
dated values  through 
the  relaxation 
procedure 


h..ix^)^l,,i>^*)['--[-^^-^-\)] 


ITN(I) 


NNP(I) 


Iteration  count  for 
the  convergence  of 
solutions  at  a  given 
NX-station. 

Number  of  points 
across  the  boundary 
layer . 


C.   HOW  TO  CHANGE  THE  ORIGINAL  PROGRAM 

Cebeci's  computer  program  is  written  entirely  in  FORTRAN 
and  consists  of  the  following  eight  FORTRAN  files: 

FILE  1    FORTRAN   Al 

FILE  2    FORTRAN   Al 

FILE  3    FORTRAN   Al 

FILE  4    FORTRAN   Al 

FILE  5 FORTRAN W- 

FILE  6    FORTRAN   Al 
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FILE  7    FORTRAN   Al 

FILE  8    FORTRAN   Al 

FILE  9    FORTRAN   Al 
Also,  the  user  needs  the  following  six  JCL  files  to  run 
changed  program: 


QUEST  1 

JCL 

Al 

ALLPDS 

JCL 

Al 

QUEST  3 

JCL 

Al 

STOPDS 

JCL 

Al 

LOADMOD 

JCL 

Al 

INTAPCLG 

JCL 

Al 

"QUEST 

1" 

lists  all 

data  se 

"ALLPDS" 


"QUEST  3" 


"STOPDS" 


System)  that  belong  to  student  user  ID 

number  and  in  group  PUB4B  (public  disk. 

volume),  which  deletes  the  data  set  180 

days  after  creation. 

allocates  space  for  a  PDS  (partitioned 

Data  Set)  on  PUB4B  prior  to  reading  a 

data  set  into  it,  using  the  utility 

program  lEFBR  14  which  pre-allocates  or 

deletes  data  sets. 

lists  the  members  of  a  PDS  and  calculates 

the  remaining  available  space  within  the 

data  set. 

places  a  FORTRAN  source  file  in  a  PDS. 
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The  SYSUT2  DD  statement  describes  the  out- 
put data  set  and  the  SYSUTl  DD  statement 
describes  the  input  data  set. 
"LOADMOD"         creates  a  load  module. 

"INTAPCLG"       compiles  50urre.  code,  iaods   ond  C^ecu/ts  iexf  file , 
A  list  of  all  JCL  files  is  given  at  the  end  of  this  chapter 
and  more  information  is  in  the  User's  Guide  to  MVS  at  NPS 
(1986). 

The  following  is  the  procedure  to  change  the  program 
and  run  it: 
Step  1 

1.  Make  changes  in  FORTRAN  files  as  required.   Remember 
those  eight  FORTRAN  files  compose  one  program. 
Thus,  the  user  must  check  all  files  which  are 
related  to  changes. 

2.  Review  and  update  FORTRAN  errors,  if  necessary. 
Step  2 

1.  Change  job  name  "xxxxxx",  user  ID  number  "nnnn"  and 
library  name  "yyy"  with  user's  own  in  all  JCL  files 
to  make  them  be  the  user's. 

2.  Submit  four  JCL  files: 

"QUESTl"   "ALLPDS"   "QUEST  3"   "LOADMOD". 

Step  3 

Do  the  following  for  all  FORTRAN  files  to  submit 

them  to  user's  library. 

1.  Open  "STOPDS"  JCL  file. 

2.  Change  the  FORTRAN  file  number  indicated  by  n  on  the 
sixth  line. 

//SYSUT2  DD  DISP  YYYLIB  (INTAIR  n) 
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3.  Insert  the  FORTRAN  file  between 

"//SYSUTl  DD,  *"  and  "/*". 

4.  Submit  "STOPDS"  JCL  file. 

5.  Check  the  result  by  any  error  message  at  the  very 
top  and  bottom  part  arid  condition  code. 

6.  Update  errors,  if  necessary. 
Step  4 

Do  the  following  to  run  the  changed  program. 

1.  Open  "INTAPCLG"  JCL  file. 

2.  Insert  the  input  data  between 

"//GO.   SYSIN   DD  ^    "  and  "/*". 

3.  Submit  "INTAPCLG"  JCL  file 

From  the  second  change,  do  steps  1,  3  and  4. 

D.   APPLICATION  OF  CEBECI ' S  PROGRAM 

Cebeci's  interactive  program  was  applied  to  a  single 
airfoil,  FX  63  -137,  at  three  Reynolds  numbers.   The  print 
flag  is  0  and  the  number  of  airfoil  coordinates  is  49  for 
all  cases.   The  results  are  compared  with  experimental  data 
and  the  turbulence  model  has  been  changed  to  get  better 
results.   The  experimental  results  are  taken  from  Reference 
9. 

Figure  4.2  and  4.3  show  the  comparison  of  lift  and  drag 
coefficients  for  three  Reynolds  numbers.   Drag  coefficients 
are  computed  at  the  trailing  edge.   It  is  seen  that  the 
results  are  closer  to  the  experimental  data  as  the  Reynolds 
number  increases. 
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Figures  4.4,  4.5  and  4.6  show  the  skin  friction 
coefficients  on  the  upper  and  lower  surfaces  as  a  function 
of  angle  of  attack  for  the  same  three  Reynolds  numbers. 
Legends  indicate  the  angle  of  attack.   As  the  Reynolds 
number  increases,  the  length  of  the  separation  bubble 
decreases  on  both  surfaces  for  the  same  angle  of  attack. 
Especially,  Figure  4.4,  which  has  a  different  label  of 
y-axis  from  the  other  figures,  indicates  an  unrealistic  flow 
on  the  upper  surface  at  low  Reynolds  number.   In  Figure  4.5 
and  4.6,  however,  the  flow  on  the  upper  surface  is  separated 
from  about  80%  chord  at  all  angles  of  attack. 

Figures  4.7,  4.8  and  4.9  present  how  the  displacement 
thickness  varies  according  to  the  angle  of  attack  and  the 
Reynolds  number.   As  the  Reynolds  number  increases,  the 
displacement  thickness  decreases  on  both  surfaces.   As  the 
angle  of  attack  increases,  it  also  increases  on  the  upper 
surface,  but  decreases  on  the  lower  surface. 

Table  4.1  represents  the  computed  and  fixed  transition 
locations  which  are  used  in  Figures  4.10   4.13.   The 
computed  transition  locations  XTRL  =  0.288922  at  ALPD  =  8 
degrees  and  XTRL  =  0.206866  at  ALPD  =  10  degrees  were 
printed  incorrectly.   In  this  case,  the  correct  transition 
locations  can  be  obtained  by  the  following  procedure: 

1.  Calculate  the  difference  of  XTRL's  at  ALPD  =  4 
degrees  and  6  degrees. 

2.  Add  the  difference  to  XTRL  of  ALPD  =  6  degrees  and 
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take  the  result  as  a  temporary  XTRL  for  ALPD  = 
8  degrees. 

3.  Fix  transition  location  with  the  temporary  XTRL 

4.  Run  the  program  to  print  XC(I),  GAMTR(I)  from 
SUBROUTINE  OUTPUT  in  FILE  3  FORTRAN  Al. 

6.   Take  XC(I),  where  the  value  of  GAMTR(I)  is  finally 
equal  to  zero,  as  the  correct  transition  location. 

6.   From  the  same  procedures  with  XTRL's  at  ALPD  =  6 

degrees  and  8  degrees,  obtain  the  correct  transition 
location  for  ALPD  =  10  degrees. 

Figures  4.10-^4.13  show  the  effect  of  variations  in  the 
empirical  constant,  Gy^        ,    in  turbulence  model  at  R^  =  0.28 
X  10   .   In  the  legend  of  Figures  4.10  and  4.11,  the  three 
numbers  are  the  values  of  the  empirical  constant,  C  means 
that  the  transition  location  is  computed  by  Eq .  (4.1)  and  F 
means  that  it  is  fixed.   The  lift  curve  is  closer  to  the 
experimental  data  as  the  empirical  constant  decreases,  but 
the  drag  curve  is  closest  when  the  empirical  constant  is 
120. 

The  reason  why  we  can  obtain  better  results  by  reducing 
the  empirical  constant  is  shown  in  Figure  4.12  and  4.13. 
GAMMATR  (  "f-i^)    is  the  intermi ttency  factor  which  is  a 
function  of  the  x-coordinate  with  values  0.0  at  the 
beginning  point  of  transition  and  1.0  at  the  ending  point  of 
transition.   The  intermittency  factor  makes  it  possible  to 
avoid  a  sudden  transition  from  laminar  to  turbulent  by 
smoothing  out  the  step-shaped  change  of  viscosity  from 
kinematic  to  eddy.   The  relation  between  the  empirical 
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constant  and  the  intermittency  factor  is  given  as: 

tr 


/^P  =      1    -   exp 


yj  tr)    Ue    R<?  ^ 


dx 


Xf^    U? 


where,  x^-  denotes  the  beginning  point  of  transition. 
The  transition  length  is  mainly  determined  by  the 
empirical  constant  Gv^   which  is  set  as  1200  in  Cebeci's 
original  program.   Decreasing  the  value  of  the  empirical 
constant  reduces  the  transition  length.   In  this  way,  the 
unrealistic  flow  shown  in  Figure  4.4  can  be  avoided  and  more 
reasonable  results,  which  are  closer  to  the  experimental 
data,  can  be  obtained.   In  Figure  4.12  and  4.13,  the  legend 
indicates  the  angle  of  attack,  and  each  of  four  line 

patterns  is  used  twice  for  two  angles  of  attack  (e.g.  ; 

-  4  degrees  and  4  degrees).   The  direction  of  curves,  as  the 
angle  of  attack  increases,  is  from  right  to  left  at  the 
upper  surface  (Figure  4.12)  and  from  left  to  right  at  the 
lower  surface   (Figure  4.13). 
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List  of  JCL  Files 


xxxxxxx  ALLPDS  JCL   Al   xxxxxxx 


//XXXXXX    JOB  (NNNN, 9999), 'ALLOCATE  PDS«,CLASS=A 

//XMAIN  0RG  =  NPGVM1  .NMNfJP 

//        EXEC   PGM=IEFBR1'^ 

//SYSPRINT  DD  SYSOUT=A 

//DDl  DD  UNIT  =  3330V,MSVGP  =  PUB'4B,DISP  =  (NEW,CATLG, DELETE), 

//        SPACE=(CYL,(^,'H,6)),DSN  =  MSS.SNNMN.YYYLIB 

// 


xxxxxxx  QUESTl   JCL   Al   xxxxKxx 


//XXXXXX    JOB  (NNNN, 9999), 'QUESTION', CLASS=A 
//XMAIN  ORG  =  NPGVMl.NrUJNP 
//   EXEC  PGM=IDCAMS 
//SYSPRINT  DD  SYSOUT=A 
//SYSIN     DD  ^ 
LISTDSET  GROUP(PUB'iB)  LEVEL  (MSS  .  SNNNN) 

/K 
// 


xxxxxxx   QUEST3   JCL   Al   xxxxxxx 


//XXXXXX    JOB  (NNNN, 9999), 'LIST', CLASS=A 

//XMAIN  0RG=NPGVM1 .NNNNP 

//   EXEC  PGM=IEHLIST 

//SYSPRINT  DD  SYSOUT=A 

//DDl  DD  UNIT=3330V,VOL=SER=MS0005,DISP=SHR 

//SYSIN     DD  X 
LISTVTOC  FORMAT, VOL = 3330 V=MSOO 05, DSNAME=MSS. SNNNN. YYYL IB 
LISTPDS   VOL =3330V=MS0  0 05, DSNAME=MSS. SNNNN. YYYL IB 

/x 

// 


xxxxxxx  STOPDS   JCL   Al   xxxxxxx 

//XXXXXX     JOB  (NNNN, 9999), 'PLACE  PDS' 

//XMAIN  0RG=NPGVM1 .NNNNP 

//       EXEC   PGM=IEBGENER 

//SYSPRINT  DD  SYSOUT=A 

//SYSIN     DD  DUMMY 

//SYSUT2    DD  DISP=(0LD, KEEP), DSNAME=MSS. SNNNN. YYYLIBdNTAIR 

//  DCB=(RECFM=FB,LRECL=80,BLKSIZE=8000) 

//SYSUTl    DD  X 

/x 

// 
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xxxxxxx   LOADMOD  JCL   Al   xxxxxxx 


//xxxxx 

JOB 

(NNNN,9999), 'CREATE  LOADMODULE' ,CLASS=C 

//KMAIN 

ORG  =  nPGVMl  .NMNfJP 

// 

EXEC 

FORTVCL,PARM.FORT='LVL(77),NOS,NOX,NOMAP' 

//FORT. 

SYSPRINT 

DD  DUMMY 

//FORT. 

SYS  IN 

DD  DISP=SHR,DSN=MSS.SHNNN.YYYLIB(INTAIR1) 

// 

DD  DISP=SHR,DSN=MSS.SNNMN.YYYLIB(INTAIR2) 

// 

DD  DISP=SHR,DSN=MSS.SMMNH.YYYLIB(INTAIR3) 

// 

DD  DISP  =  SHR,DSN  =  MSS.SNriNfJ.YYYLIB(INTAIR^) 

// 

DD  DISP=SHR,DSN=MSS.SNfU<N.YYYLIB(INTAIR6) 

// 

DD  DISP  =  SHR,DSN  =  M3S.SNNfiN.YYYLIB(INTAIR7) 

// 

DD  DISP=SHR,DSN=MSS.SNNfUi.YYYLIB(INTAIR8) 

// 

DD  DISP  =  SHR,DSN  =  MSS.SNfJNN.YYYLIB(INTAIR9) 

//LKED. 

SYSLMOD 

DD  DISP=SHR,DSN=MSS.SNNNN.YYYLIB(INTCAS) 

//LKED. 
// 

SYSIN  DC 

DSN=NULLFILE 

xxxxxxx   INTAPCLG  JCL   Al   xxxxxxx 


//XXXXXX    J 
//XMAIN  ORG=N 
//        EXEC 
//        REGI 
//FORT.SYSPRI 
//FORT. SYSIN 
// 
// 
// 
// 
// 
// 
// 

//LKED.SYSPRI 
//LKED. SYSIN 
//GO.FTOIFOOI 
//GO.FT02F001 
FT03F001 
FTO^FOOl 
FT99F001 
FT08F001 
FT09F001 
FTlOFOOl 
FTUFOOl 
FT06F001 
SYSIN 


//GO, 

//GO, 

//GO, 

//GO, 

//GO 

//GO 

//GO 

//GO 

//GO 

/x 

// 


OB  (NNNN,9999), 
PGVMl .NNNNP 

FORTVCLG.PARM 
ON.GO=1280K 
NT  DD  DUMMY 

DD  DISP=SHR, 
DD  DISP=SHR, 
DD  DISP=SHR, 
DD  DISP=SHR, 
DD  DISP=SHR, 
DD  DISP=SHR, 
DD  DISP=SHR, 
DD  DISP=SHR, 
NT  DD  DUMMY 

DD  DSN=NULLF 
DD  UNIT=SYSDA, 
DD  UNIT=SYSDA, 
DD  UNIT=SYSDA, 
DD  UNIT=SYSDA, 
DD  UNIT=SYSDA, 
DD  UNIT=SYSDA, 
DD  UNIT=SYSDA, 
DD  UNIT=SYSDA, 
DD  UNIT=SYSDA, 
DD  SYSOUT=A 
DD  X 


'COMP  LOAD  AND  GO* 
.F0RT='LVL(77),N0S 


DSN  = 

DSN  = 
DSN  = 
DSN  = 
DSN  = 
DSN  = 
DSN  = 
DSN  = 


MSS 
MSS 
MSS 
MSS 
MSS 
MSS 
MSS 
MSS 


SNNNN 

SNNfJN 

SNNNN 

SNNNN 

,  SNNNfJ 

SNNNN 

,  SNNfUJ 

,  SNNNN 


.YYYL 
.YYYL 
.YYYL 
.YYYL 
.YYYL 
.YYYL 
.YYYL 
.YYYL 


,CLASS=G 

, N0X,N0MAP', 


IB(INTAIRl) 
IBdNTAIRZ) 
IB(INTAIR3) 
IB(INTAIR^) 
IB(INTAIR6) 
IB(INTAIR7) 
IBCINTAIRS) 
IB(INTAIR9) 


ILE 

SPACE=(CYL,(1,1 

SPACE=(CYL,(1,1 

SPACE=(CYL,(1,1 

SPACE=(CYL,(1,1 

SPACE=(CYL,(1,1 

SPACE=(CYL,(1,1 

SPACE=(CYL,(1,1 

SPACE=(CYL,(1,1 

SPACE=(CYL,(1,1 
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TABLE  4.1    TRANSITION  LOCATION 


COMPUTED 

FIXED 

ALPD 

XTRL 

XTRU 

XTRL 

XTRU 

-4 

0.010680 

0.584178 

0.010680 

0.584178 

-2 

0.326179 

0.542882 

0.326179 

0.542882 

0 

0.391919 

0.501232 

0.391919 

0.501232 

2 

0.458252 

0.478890 

0.458252 

0.459243 

4 

0.524927 

0.416490 

0.524927 

0.416490 

6 

0.652658 

0.353115 

0.632385 

0.353115 

8 

0.288922 

0.288922 

0.676143 

0.288922 

10 

0.206866 

0.224747 

0.718631 

0.206866 

80 


C\2 
CO 

CM 

yl' 

■/// 

/ 

-^ 

CO 

J 

/ 

/ 

/ 

d 

/^ 

// 

/  / 

r  / 

/ 

-* 

LECxEND 
EXPERIMENTAL 

d 

// 

/ 

1200  ( 

1 

20  F 

o 

__12  F_ 

d  J 

1          I          1          1 

-6.0  -4.0  -2.0   0.0   2.0   4.0   6.0   8.0   10.0  12.0 


ALPHA  (DECxREE) 

Figure  4.10   Comparison  Of  Lift  Curves 

For  Different  Empirical 

Constants  G^     (  Re  =  0.28  x  10   ) 
"tr 


81 


CO 

C\] 

^ 

'/^       "^ 

^       ^   ^ 

X 

■"^  ^ - -' ' 

1—1 
CO 

// 

/// 
//  ''  ' 

/ 
/ 
/ 

/ 

d 

( 

' '  '  / 

1 

\     \ 

\    ^ 
vN^ 

^ 

I  Fr.Fivn 

d 

X\ 

\^- 

EXPERIMENTAL 

\ 

1200  C 

120  F 

o 

12  F 



d  J 

— ^ ^ ^ 1 

0.008    0.010         0.012         0.014         0.016         0.018         0.020 

CD 


Figure  4.11   Comparison  Of  Drag  Curves 
For  Different  Empirical 


Constants  G 


7 

*  tr 


(  Re  =  0.28  X  10   ) 


82 


CJ 


LEGEND 
-4 

w 

.0 

0 
"2"" 
4 

6 
8 

"lo' 

0.0                0 

I — '■ — '—  "^I                 1                 1 

2                  0.4                  0.6                  0.0                  I 

x/c 

< 


(a)   G, 


1200 


tr 


0.0       0.2       0.4       0.6 

X/C 

(b)   G     =  120 
•tr 


0.8 


<        ■     1 

/  //, 

■■// 

/   / 

o 

/       /    / 

// 

1.0 


/     :'    / 

/  ■■'/ 

:  \  \ 

\'4 

\  . ' 

vA 

0.0 

M 

\\\\ 

0.0      0.2       0.4       0.6       0.8       1.0 

X/C 

(c)   G     =  12 
«  tr 


Figure  4.12    Comparison  Of  Transition  Length 

For  Different  Empirical  Coefficients  G.^ 

6    '^tr 
At  Upper  Surface   (  Re  =  0.28  x  10   ) 


83 


< 


< 


-M  ~ 

y' 

, 

X              / 

y             ^ 

^ 

^ 

/ 

y 

X 

,/                ' 

/ 

/ 

^/ 

X 

o 

/ 

y  y 

^.^ 

o 

i 

< 


u 


0.0 


0.0      0.2       0.4       0.6 

x/c 

(a)   G  .   =  1200 
Ttr 


0.8 


0.0      0.2       0.4       0.6 

X/C 

(b)   G     =  120 
4  tr 


0  8 


0.2 


0.4 


0.6 


0.8 


X/C 


1.0 


1.0 


/  :  /    ;     / 

/■/.. 

v'\\\ 

/// 

/   ;;    1 

W^ 

o 

/ '  /  / 

//^ 

1.0 


LEGEND 

-4 

-2 

0 
"Z 

4 


LEGEND 
-4 


—  O 


0 
2 
4 


_6_ 
8 

io 


LEGEND 
-4 


_o 


_6_ 
0 

io 


(c)   G^,   =  12 


Figure  4.13    Comparison  Of  Transition  Length 

For  Different  Empirical  Coefficients  G  . 

6   'tr 
At  Lower  Surface   (  Re  =  0.28  x  10   ) 


8A 


V.   CONCLUSIONS 

Cebeci's  interactive  computer  program  was  applied  to 
the  Wor tmann-Al thaus  FX  63-137  airfoil  to  show  the 
capability  of  strong  v iscous/ inviscid  interaction  methods  to 
predict  airfoil  flows  at  low  Reynolds  numbers. 

From  the  comparisons  with  the  experimental  data,  it 
was  confirmed  that  the  results  are  closer  to  the 
experimental  data  as  the  Reynolds  number  increases. 

Also,  much  better  results  were  obtained  by  decreasing 
the  empirical  constant  G^ 

Therefore,  it  was  concluded  that  the  boundary  layer 
transition  model  has  an  important  influence  on  the 
predictive  capability  of  viscous/ inviscid  interaction 
methods. 
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